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Abstract 

In  the  present  paper,  a  three-dimensional,  time-dependent  SOFC  numerical  model  is  defined,  considering  all  the  phenomena  occurring  in  each 
component  of  the  fuel  cell.  All  the  equations  are  written  in  a  partial  differential  form,  thus  the  model  is  independent  from  the  cell  geometry  (planar, 
tubular,  monolithic)  and  the  modeling  approach  (i.e.  time-dependent,  3D,  2D).  The  triple-phase-boundary  (TPB)  is  modeled  as  a  finite  region,  as 
well  as  a  non-dimensional  interface  area.  In  the  latter  case  the  source  terms  of  the  energy,  mass  and  electrical  conservation  equations  are  treated  as 
a  boundary  conditions.  Finally,  a  literature  review  is  conducted.  Some  results  from  the  models  previously  developed  by  the  authors  or  found  in  the 
literature  survey  are  presented.  Each  of  them  can  be  reviewed  as  a  simplified  version  of  the  model  presented  in  the  present  paper.  The  simplified 
assumptions  and  the  related  omissions  that  each  of  those  approaches  imply  are  also  analyzed. 

©  2005  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

Solid  oxide  fuel  cells  (SOFCs)  are  considered  promising 
energy  conversion  devices  thanks  to  their  several  potential  ben¬ 
efits,  including  low  pollutant  emissions,  high-energy  efficiency, 
the  possibility  of  using  different  kinds  of  fuels  and  the  possibility 
to  build  CHP  and  hybrid  systems. 

However,  SOFC  technology  is  still  in  embryonic  infancy 
and  many  problems  (i.e.  mechanical  stress,  electrode  sin¬ 
tering,  electrode  and  interconnect  materials  and  fabrication, 
startup  time)  must  be  solved  in  order  to  achieve  the  goal  of 
a  highly  efficient  and  clean  energy  system  with  at  least  the 
same  reliability,  costs  and  lifetime  of  the  ‘Traditional”  energy 
systems. 

This  challenge  entails  an  increased  research  effort  for  indus¬ 
try  and  research  institutions  in  order  to  develop  advanced  numer¬ 
ical  models  and  computer  codes,  that  can  be  used  in  both 
industrial  and  research  environments  for  fuel  cells  design  and 
development.  This  could  yield  a  better  understanding  of  the 
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physical  processes  and  help  fuel  cell  designers  to  define  more 
promising  strategies. 

The  performance,  the  reliability  and  the  duration  of  an  SOFC 
is  directly  related  to  the  processes  that  occur  within  the  cell.  The 
control  of  these  processes  requires  primarily  an  understanding 
of  SOFC  physics  and  chemistry  as  well  as  the  capability  to  mod¬ 
ify  one  or  more  interdependent  processes  parameters  in  a  given 
direction.  This  task  is  made  even  more  difficult  by  the  complex¬ 
ity  and  the  vast  number  of  physical  and  chemical  processes  that 
occur  simultaneously  in  an  SOFC  and/or  in  an  energy  system 
based  on  SOFCs.  All  this  along  with  the  high  number  of  param¬ 
eters  involved  necessarily  make  computer  models  attractive  for 
their  potential  of  modeling  physical  fluid  phenomena  that  cannot 
be  easily  measured. 

Enormous  progress  has  been  made  in  the  last  decade  on 
numerical,  analytical  and  computational  tools  for  SOFCs.  One  of 
the  first  significant  publications  on  SOFC  modeling  is  Bossel’s 
[1],  published  in  1992  as  a  result  of  an  Electric  Power  Research 
Institute  (EPRI)  project.  The  analysis  assesses  the  main  differ¬ 
ences,  in  terms  of  performance,  among  different  SOFC  config¬ 
urations,  using  the  same  materials  characteristic.  The  model  is 
quite  simple,  but  at  the  same  time  provides  a  very  important  tool 
to  investigate  cross-plane  and  in-plane  resistance.  The  SOFC 
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configurations  studied  at  that  time  still  constitutes  the  main  geo¬ 
metrical  design  and  the  results  of  that  study  still  represent  a 
relevant  base  point  for  understanding  the  performance  potential 
of  any  SOFC  configuration.  The  model  developed  by  Achenbach 
[2],  instead,  is  based  on  a  three-dimensional  and  time-dependent 
approach  for  planar  solid  oxide  fuel  cell  simulation.  This  model 
appears  to  be  quite  accurate,  and  uses  differential  and  finite 
equations  that  allow  detecting  gas  concentration,  current  den¬ 
sity,  temperature  variation  within  the  cell,  along  with  the  effects 
of  different  fuels,  flow  manifolding  and  of  three-dimensional 
phenomena.  The  model  considers  the  effect  of  ohmic  losses 
and  anode  and  cathode  polarizations,  in  the  form  of  “activation” 
polarization. 

The  approach  used  by  Bessette  et  al.  [3]  for  a  tubular  solid 
oxide  fuel  cell  model,  is  similar  to  Achenbach’s.  Besides  the 
ohmic  and  activation  losses,  a  method  to  evaluate  concentration 
polarization  is  given.  The  results  show  the  fuel  cell  local  char¬ 
acteristics  for  different  axial  and  angular  positions.  Moreover, 
the  effects  of  different  oxidants  and  fuels  are  shown. 

Several  other  publications  can  be  found  in  the  literature, 
which  illustrate  more  accurate  models,  thanks  to  the  enhanced 
computational  performance.  Examples  of  very  detailed  models 
are  provided  by  Ferguson  [4]  or  Petruzzi  et  al.  [5]. 

However,  considerable  work  is  still  necessary  in  developing 
and  testing  sub-models  to  describe  unresolved  individual  physi¬ 
cal  processes  that,  especially  for  the  porous  electrodes,  introduce 
an  excess  of  empiricism  into  computation. 

Numerical  modeling  is  now  playing  an  important  role  and 
represents  a  critical  aspect  of  the  technology  development  pro¬ 
cess  in  industry.  This  is  mainly  due  to  the  increase  in  competition 
in  all  fields  of  engineering  and  to  the  enormous  improvements  in 
computer  technology  in  the  last  two  decades.  Moreover,  model¬ 
ing  is  widely  utilized  within  industry  because  it  is  a  valid  tool  to 
investigate  physical  fluid  systems  more  cost  effectively  and  more 
rapidly  than  with  experimental  procedures.  In  particular,  the 
process  of  modeling  is  integrated  into  the  entire  engineering  pro¬ 
cess  that  usually  foresees  a  combined  experimental-numerical 
approach: 

-  numerical  simulations  are  necessary  to  have  direct  “insight” 
into  the  system  and  to  allow  studying  different  components 
design,  stack  and  cell  configurations  and  layouts,  fuels,  etc. 
and  to  determine  optimum  operating  conditions;  this  reduces 
the  number  of  experimental  tests  to  be  performed  and  the 
costly  and  time  consuming  physical  prototyping;  in  brief,  three 
roles  can  be  identified  for  numerical  modeling:  interpretation, 
design  and  prediction; 

-  measurements  are  necessary  to  validate  numerical  models 
and  to  make  the  final  choice  among  those  few  configurations 
selected  through  computational  simulations. 

However,  numerical  modeling  must  be  used  very  carefully, 
especially  when  it  is  used  for  prediction  purposes.  It  must  be 
borne  in  mind  that  numerical  modeling  means  reproducing  the 
behavior  of  a  system  by  solving  a  set  of  equations  describing  the 
evolution  of  variables  that  define  the  state  of  the  system.  This 
evolution  is  governed  by  processes  that  are  extremely  complex 


and  some  of  them  are  completely  unknown.  This  means  that 
even  with  a  high  level  of  model  complexity,  models  are  only  a 
simplified  representation  of  real  physics.  Therefore,  the  accurate 
prediction  of  actual  performances  cannot  be  guaranteed  and  is 
instead  very  difficult  even  after  appropriate  validation  of  the 
models. 

The  present  paper  is  organized  in  order  to  follow  the  usual 
steps  in  the  development  of  a  numerical  model: 

-  understanding  the  physical  system  and  translating  it  into  math¬ 
ematical  equations;  in  their  most  general  form,  these  are 
usually  a  set  of  partial  differential  equations  (mathematical 
model); 

-  these  equations  cannot  be  usually  analytically  solved  and 
are  discretized  to  allow  numerical  solution;  spatial  domain 
is  divided  into  small  elements  constituting  the  grid  or  mesh 
(numerical  model); 

-  after  generating  a  reasonably  fine  numerical  grid  on  which 
the  discrete  algebraic  equations  will  be  solved,  a  set  of 
problem-dependent  boundary  and  initial  conditions  are  spec¬ 
ified  (boundary  conditions); 

-  it  is  then  shown  how  the  complete,  3D  and  time-dependent 
problem  can  be  simplified,  according  to  the  specific  needs,  to 
steady  state  conditions,  2D,  ID  and  OD; 

-  once  a  computational  model  has  been  developed  for  a  partic¬ 
ular  application,  the  results  must  be  compared  with  physical 
reality;  this  validates  the  code  and  avoids  blind  acceptance  of 
numerical  results,  which  is  often  a  problem  with  commercial 
software  use  (validation); 

-  some  results  obtained  by  the  authors  in  previous  studies, 
as  well  as  simplified  applications  of  the  model  are  shown 
(results). 

2.  Mathematical  model 

A  comprehensive  analysis  of  solid  oxide  fuel  cells  (SOFC) 
phenomena  requires  an  effective  multidisciplinary  approach. 
Chemical  reactions,  electrical  conduction,  ionic  conduction  and 
heat  transfer  take  place  all  at  the  same  time  and  are  tightly  cou¬ 
pled. 

SOFC  can  be  manufactured  in  different  geometrical  config¬ 
urations,  i.e.  planar,  tubular  or  monolithic. 

Regardless  of  the  geometrical  configuration,  a  solid  oxide 
fuel  cell  is  always  composed  of  two  porous  electrodes  (anode 
and  cathode),  a  dense  electrolyte,  an  anodic  and  a  cathodic  gas 
channel  and  two  current  collectors.  For  sake  of  simplicity,  in 
Fig.  1,  the  planar  configuration  is  taken  as  reference.  It  is  inter¬ 
esting  to  note  that  for  the  flat  planar  configuration  (Fig.  1),  the 
current  collector  acts  also  as  a  gas  channel,  while  for  the  tubular 
and  monolithic,  one  electrode  acts  as  the  gas  distributor. 

The  task  of  the  gas  channel  is  to  guide  the  gas  on  the  electrode 
surface  for  enabling  a  uniform  distribution.  If  the  cell  operates  on 
a  hydrocarbon,  rather  than  on  hydrogen,  the  reforming  reaction 
(1)  takes  place  in  the  anodic  gas  channel: 

CnUmOp  +  in  -  p) H2O  =  nC O  +  (n  -  p  +  m/2)H2 


(1) 
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Fig.  1.  Schematic  representation  of  a  flat  planar  SOFC. 


•  Chemical  reactions  (1)  and  (2),  taking  place  in  the  anodic  gas 
channel  if  hydrocarbons  are  internally  reformed.  The  steam 
reforming  reaction  (1)  is  endothermic,  thus  heat  is  absorbed 
from  the  surroundings. 

•  Ohmic  resistance;  the  ionic  and  electronic  currents,  when 
passing  through  the  electrolyte  and  the  electrodes,  generate 
the  so-called  Joule  heating. 

•  Activation  over-potential  due  to  the  activation  energy;  the  cur¬ 
rent  flow  generates  an  electrochemical  loss  that  is  translated 
into  a  heat  source. 


Furthermore,  carbon  monoxide  reacts  with  water  to  generate 
additional  hydrogen  and  carbon  dioxide: 

CO  +  H20  C02  +  H2  (2) 

Since  the  anode  is  porous,  the  gases  in  the  anodic  gas  channel 
permeate  it.  The  hydrogen  generated  from  reactions  (1)  and  (2) 
reacts  with  the  oxygen  ions  coming  from  the  cathode  through 
the  electrolyte,  thus  generating  water: 

H2  +  O2-  -*  H20  +  2e~  (3) 

The  water  produced  goes  back  to  the  anode  channel  and  exits 
the  cell,  together  with  the  other  unoxidized  gas  species.  The 
result  for  the  anodic  side  is  a  hydrogen  flux  from  the  channel  to 
the  anode-electrolyte  interface  and  a  water  flux  on  the  opposite 
path. 

At  the  cathode  side,  instead,  oxygen  in  the  cathodic  channel 
permeates  the  cathode  where  ions  form  and  migrate  through  the 
electrolyte: 

i02  +  2e-  O2-  (4) 

The  result  for  the  cathode  is  a  one  way  flux  from  the  channel 
to  the  electrode-electrolyte  interface  (Fig.  2). 

The  electrons  released  from  Eq.  (3)  flow  from  the  reaction 
zone  (typically  very  close  to  the  anode/electrolyte  interface)  to 
the  current  collector,  through  the  anode.  At  the  cathode  side, 
these  electrons  flow  from  the  current  collector  to  the  cathodic 
reaction  zone,  where  reaction  (4)  takes  place. 

Heat  is  generated  inside  the  cell  by  the  following  heat  sources: 

•  Electrochemical  reactions  (3)  and  (4),  taking  place  very 
close  to  the  electrode-electrolyte  interface  (i.e.  at  the  so- 
called  triple-phase-boundary  (TPB))  generate  heat,  due  to  the 
entropy  increase. 


The  heat  generated  inside  the  cell  is  removed  by  the  incoming 
anodic  and  cathodic  gases.  The  heat  transfer  between  the  cell 
components,  the  gases  and  the  surroundings,  acts  for  conduction, 
convection  and  radiation. 

In  the  following,  the  mathematical  model  of  each  component 
is  defined,  thus  obtaining  the  system  of  equations  that  analyti¬ 
cally  describes  a  generic  SOFC.  Furthermore,  the  completion  of 
the  model  with  the  boundary  conditions  is  discussed. 

2.7.  Channels  flow 

Using  an  Eulerian  approach,  the  description  of  fluid  motion 
requires  the  determination  of  the  thermodynamic  state,  in  terms 
of  sensible  fluid  properties,  pressure,  P ,  density,  p  and  temper¬ 
ature,  T,  and  of  the  velocity  field  u  ( x ,  t)  [6-10]. 

Therefore,  in  a  three-dimensional  space  for  a  given  thermo¬ 
dynamic  system  having  two  intensive  degrees  of  freedom,  we 
have  six  independent  variables  as  unknowns  of  the  so-called 
“thermo-fluid  dynamic  problem”,  thus  requiring  six  indepen¬ 
dent  equations.  The  six  equations  are  given  by  the  equation  of 
state  and  the  three  fundamental  principles  of  conservation: 

-  equation  of  state  for  involved  gas  species; 

-  mass  continuity;  actually,  if  a  mixture  of  gases  is  present, 
species  mass  fractions  are  additional  unknowns,  thus  requiring 
the  mass  conservation  equation  for  each  specie; 

-  Newton  second  law  or  momentum  equation  (three  equations 
in  a  three-dimensional  space  x,  y  and  z); 

-  energy  conservation  (First  law  of  thermodynamics). 

Since  gases  in  SOFCs  are  far  from  the  critic  conditions,  per¬ 
fect  gas  equation  of  state  is  usually  employed: 

P  =  pRT  (5) 


H2,  CO.  CO2, 


H20.  H2,  CO. 


N2,  etc 


O2,  N2 


C02;_N2,  etc 


02,  N2 

■» 


Fig.  2.  Gas  fluxes  inside  the  SOFC. 
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where  P  is  the  gas  pressure,  p  the  density,  R  the  universal  gas 
constant  and  T  is  the  temperature. 

The  conservation  equations  can  be  expressed  in  either  a  dif¬ 
ferential  or  integral  form.  The  differential  form,  usually  called 
primitive  formulation,  is  the  most  used  in  fluid  mechanics,  and 
is  used  in  the  present  formulation. 

Considering  that  the  working  fluid  in  the  gas  channels  of  an 
SOFC  is  a  reacting  mixture  of  gases,  mass  conservation  can  be 
analytically  described  by  the  following  transportation  equation: 

dpYi  _  - 

— - h  u  •  V(pYi)  =  Vm;  +  (Oi  (6) 

ot 

where  i  denotes  the  generic  ith  specie,  coi  the  rate  of  production 
or  consumption  of  the  specie,  mz-  the  mass  diffusion  flux,  7/  the 
mass  fraction  of  the  ith  species,  and  u  is  the  gas  velocity.  Accord¬ 
ing  to  Fick’s  law  of  binary  diffusion,  the  ith  specie  diffusion  flux 
can  be  read  as  follows: 


mi  = 


pDiVW) 


where  Df  is  the  mass  diffusion  coefficient. 

The  overall  fluid  density  is  the  sum  p  =  J2?=\YiP  over  all 
species  and  leads  to  the  standard  continuity  equation: 

dp 

—  +  V  (pu)  =  oj  (8) 

ot 

where  co  is  a  source  term  that  takes  into  account  mass  variation 
due  to  the  electrochemical  reaction. 

In  order  to  calculate  the  physical  properties  of  the  multicom¬ 
ponent  mixture,  including  viscosity,  specific  heats  at  constant 
volume  and  at  constant  pressure  and  laminar  thermal  conductiv¬ 
ity,  the  assumption  that  the  components  form  an  ideal  mixture 
is  usually  made.  An  arbitrary  constitutive  fluid  property  may  be 
calculated  from  the  property  value  for  the  fluid  ith  component 
through  the  following  equation: 

J  =  T>  Yj  (9) 


cosity  and  thermal  conductivity  on  the  fluid  state.  However,  in 
many  situations  they  are  assumed  to  be  constant. 

In  the  above  formulation  of  the  governing  equations,  the 
velocity  field  and  the  pressure  field  are  treated  as  unknowns 
in  the  fluid  domain. 

It  is  important  to  notice  that  the  general  equations  of  fluid 
motion  are  independent  of  the  coordinate  system.  Therefore,  Eq. 
(10),  given  in  the  Cartesian  coordinate  system,  can  be  converted 
and  used  with  a  specific  coordinate  system  orientation  relative 
to  the  flow  field  (i.e.  tubular  SOFC). 

If  the  divergence  of  Eq.  (10)  is  taken,  the  pressure  field  can 
be  calculated  as  a  solution  of  Poisson’s  equation: 

V2P  = -V(mVm).  (11) 

The  physics  of  laminar  and  turbulent  incompressible 
flows  are  well  represented  by  the  above  formulation  of  the 
Navier-Stokes  equations.  Closure  of  this  system  of  equations 
is  a  non-trivial  task  and  closed  solutions  are  available  only  for 
exceptional  cases.  Considering  that  the  gas  speed  in  SOFC  gas 
channels  is  always  very  low,  it  is  a  common  practice  to  assume 
a  laminar  flow  [11]. 

This  assumption  significantly  reduces  the  computational  cost 
and  allows  to  quickly  determine  the  velocity  profile  and  the 
pressure  drop  through  the  channel,  AP  oc  pu^/ 2,  being  um  the 
average  channel  velocity. 

Temperature  profiles  inside  the  cell  can  be  obtained  through 
the  energy  equation,  which  can  be  written  in  numerous  but  equiv¬ 
alent  forms.  Indicating  with  e  the  fluid  energy  per  unit  mass,  the 
following  equation  for  the  conservation  of  energy  is  derived: 

de 

p—  +pu-'Ve  =  VQ  +  Sq  (12) 

dt 

where  Sq  is  the  volumetric  heat  source  term  and  Q  is  the  heat 
flux  vector  only  from  conduction,  which,  according  to  Fourier’s 
law,  can  be  related  to  the  temperature  field: 

Q  =  -AVT  (13) 


Regarding  the  momentum  conservation,  the  Navier-Stokes 
equations,  governing  unsteady  compressible  viscous  flows,  are 
the  basis  for  all  models.  However,  the  thermal  and  chemically 
reacting  flows  through  SOFC  gas  channels  are  often  charac¬ 
terized  by  relatively  slow  flow  motion,  much  slower  than  the 
speed  of  sound,  and  by  partially  varying  mixture  density.  Varia¬ 
tions  of  the  density  are  caused  in  these  cases  either  internally  by 
heat  release  of  chemical  reactions  or  externally  by  wall  heating 
and  by  mass  variation  due  to  the  electrochemical  reaction  but 
not  by  compression  of  the  gas.  This  situation  is  characterized 
by  a  low  Mach  number,  (near  the  incompressible  limit)  and  is 
consequently  modeled  by  a  scale  analysis  of  the  Navier-Stokes 
momentum  equations  with  respect  to  small  Mach  numbers: 

- \-uWu  =  —  VP -l - V2w  +  /  (10) 

dt  p 

where  /  is  a  specific  prescribed  body  force  (i.e.  gravity  force) 
and  /x  is  the  fluid  dynamic  viscosity.  The  system  of  equations 
should  still  be  supplemented  by  the  laws  of  dependence  of  vis¬ 


where  A  is  the  thermal  conductivity  coefficient. 

The  fluid  energy  per  unit  mass,  in  its  most  general  form,  is 
the  sum  of  the  specific  internal  energy  and  the  specific  kinetic 
energy: 

u2  _  _> 

e  =  CvT+  —  +g-r  (14) 

where  cv  is  the  specific  heat  at  constant  volume,  g  the  gravity 
acceleration  and  r  is  the  unity  gravity  vector. 

Considering  that  the  working  fluid  is  a  mixture  of  gases,  the 
potential  energy,  g  •  r  and  the  kinetic  energy,  u2/ 2,  are  usually 
neglected  and  the  fluid  energy  per  unit  mass  coincides  with  the 
specific  internal  energy. 

The  heat  source  term  Sq  in  Eq.  (12)  incorporates  the  viscous 
dissipation  (always  positive)  and  the  heat  generated  or  absorbed 
by  the  chemical  reactions  (i.e.  shift  and  reforming  reactions). 
Furthermore,  the  convective  thermal  flux  at  the  boundaries  and 
the  radiative  heat  transfer  with  solid  cell  components  must  be 
considered  (boundary  or  interface  source  terms). 
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In  any  multicomponent  system  involving  chemical  reactions, 
mass  and  energy  source  terms  must  be  considered  in  order  to 
solve  the  thermo-fluid  dynamics  problem.  These  reactions  can 
be  written  in  the  following  form: 


current  conservation: 


dpQ 

dt 


+  v-  j  = 


j  in  the  TPB 
0  elsewhere 


(17) 


Ns  Ns 

X) rikCi  O  k=  l,Nr  (15) 

i= 1  i= 1 

where  Nr  is  the  number  of  reactions,  Ns  the  number  of  species, 
rik  and  pik  the  stoichiometric  coefficient  of  the  reactants  and 
the  products  of  the  Mi  reaction  and  q  is  the  chemical  of  the 
ith  specie.  Knowing  the  reaction  rates  of  the  forward  and  the 
backward  reactions  (i.e.  Arrhenius  equation),  the  source  term 
cot  in  the  species  transport  Eqs.  (6)  can  be  calculated.  Moreover, 
each  reaction  is  characterized  by  an  absorption  (endothermic) 
or  a  release  (exothermic)  of  heat,  which  allows  to  determinate 
the  chemical  source  term  in  the  energy  equation. 

Regarding  the  gas  channels  of  an  SOFC,  significant  chemi¬ 
cal  reactions  may  take  place  only  in  the  anode  gas  channel:  the 
reforming  reaction  (1)  and  the  shift  reaction  (2).  The  shift  reac¬ 
tion  is  a  fast  process  and  is  usually  assumed  to  be  in  equilibrium. 
The  steam  reforming  process  instead  is  much  slower,  thus  the 
kinetic  of  the  reaction  needs  to  be  considered. 

2.2.  Electrodes 

To  solve  the  electrical  problem,  two  variables  need  to  be 
found,  i.e.  the  electrical  potential  (scalar)  and  the  current  den¬ 
sity  (vector).  Fig.  3  represents  the  anode  schematically.  SOFC 
electrodes  are  typically  made  of  mixed  electronic  and  ionic  con¬ 
ducting  materials  [12],  thus  both  electrons  and  ions  can  flow 
through  them.  Regarding  the  anode  (Fig.  3),  ions  coming  from 
the  electrolyte  flow  through  the  anode,  thus  combining  with  H2 
and  releasing  electrons  at  the  so-called  triple-phase-boundary 
(TPB)  zone.  Consequently,  electrons,  flow  towards  the  anodic 
current  collector.  The  electrical  problem  is  solved  with  one  vec¬ 
tor  and  one  scalar  relationship. 

The  vector  relationship  is  given  by  Ohm’s  law: 


J  =  -crV0  (16) 

where  J,  o  and  0  are,  respectively,  current  density,  conductivity 
and  electrical  potential,  and  can  be  referred  to  either  the  ionic 
or  the  electronic  current.  The  scalar  relation  is  provided  by  the 


H2O 


Tocurrent 

collector 


Triple  Phase 
Boundary 


POROUS 

ANODE 


Fig.  3.  Triple-phase-boundary  (TPB)  and  anodic  reaction. 


where  pQ  is  the  electronic  or  ionic  charge  density,  and  j  is  the 
current  generation  rate,  produced  at  the  three-phase-boundary. 
Since  the  cell  charging  time  is  very  short,  compared  to  the  other 
phenomena  acting  in  the  fuel  cell  [13],  this  can  be  considered  at 
steady  state,  i.e.  the  first  term  on  the  left  hand  side  of  Eq.  (17) 
can  be  neglected. 

The  current  source  is  regulated  by  the  Bultler-Volmer  equa¬ 
tion: 


exp 


act\ 

V  RT  ) 


—  exp 


/Q^actV 
V  RT  J. 


(18) 


In  expression  (18)  jo  is  the  exchange  current  density  and 
depends  on  the  electrode  properties  and  operation,  a\  and  a 2 
are  two  transfer  coefficients,  and  r\ act  represents  the  latter.  The 
activation  loss  can  quantified  as: 


7?  act  —  Frev  0electrode  ^electrolyte 


(19) 


where  Vrev  is  the  reversible  potential  between  the  electrode  and 
the  electrolyte,  and  can  be  computed  with  Nemst’s  equation. 

Since  the  TPB  boundary  is  usually  very  small,  compared  to 
the  electrode  and  electrolyte  thickness,  the  current  generation  is 
sometimes  considered  to  take  place  at  the  electrode-electrolyte 
boundary  [4,14,15],  thus  Eq.  (18)  is  considered  a  boundary  con¬ 
dition  for  Eq.  (17)  with  no  source  term,  and  the  term  jo  takes  into 
account  that  only  a  part  of  the  interface  is  active  in  the  reaction, 
as  explained  in  Section  4. 

Combining  Eqs.  (16)  and  (17),  the  following  relationship  is 
obtained: 


—  —  in  the  TPB 
o 

0  elsewhere 


(20) 


It  is  interesting  to  observe  that,  due  to  the  solid  state  of  the 
entire  cell,  and  due  to  the  non-perfect  planarity  of  the  com¬ 
ponents,  the  so-called  contact  resistance  is  generated  at  the 
electrode/electrolyte  interface  and  at  the  electrode/current  col¬ 
lector  interface.  Although  it  is  difficult  to  yield  a  theoretical 
formulation  for  this  loss,  experimental  or  empirical  data  can  be 
given  as  an  additional  boundary  condition. 

Considering  the  porous  media  as  continuum,  the  energy  equa¬ 
tion  can  be  written  as: 

3T 

c—  +  VQ  +  Sq=  0  (21) 

ot 

where  T  and  c  are  the  temperature  and  the  specific  heat  of  the 
medium,  respectively,  and  Q  has  already  been  defined  (Eq.  (13)). 
The  source  term  represents  the  heat  generated  due  to: 


1 .  Joule  effect,  due  to  the  ohmic  resistance. 

2.  Reversible  heat,  associated  to  chemical  reactions  (3)  or  (4). 

3.  Activation  loss. 
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While  the  Joule  effect  takes  place  in  all  the  electrode  domain, 
the  last  two  heat  sources  are  localized  in  the  TPB. 

{crV0  •  V</>  H - TAS  +  Jr/act  in  the  TPB 

J  J  2  p  /act  (22) 

<jV0  •  V0  elsewhere 

In  expression  (22),  AS  is  the  entropy  change  associated  with 
either  reaction  (3)  or  reaction  (4). 

Similarly  to  Eqs.  (17)  and  (18),  if  the  reaction  is  considered  at 
the  electrode/electrolyte  interface,  Eq.  (22)  becomes  a  boundary 
condition: 

Sq  =  aV0  •  V0  (23) 


The  diffusion  coefficients  (A',eff)  are  modified  to  take  into 
account  the  interaction  with  the  pore  walls.  The  Bruggemann 
correction  allows  the  evaluation  of  these  coefficients,  through 
the  following  expression  [17]: 

A‘,eff  =  sTDi  (27) 

where  r  is  the  tortuosity  of  the  porous  medium,  Di  the  ordinary 
diffusion  coefficient  and  A',eff  is  the  diffusion  coefficient  in  the 
porous  medium.  Some  authors,  e.g.  [20,21,27],  instead,  use  the 
following  expression: 

A,eff  =  -Di  (28) 

r 


The  heat  generation  is  considered  as  a  boundary  condition  at 
the  electrode/electrolyte  interface. 

Due  to  the  structure  of  the  electrodes,  the  velocity,  pressure 
and  species  distribution  needs  to  be  studied  using  the  theory  of 
mass  transport  in  porous  media.  This  topic  has  been  extensively 
studied  over  the  years,  and  reported  in  the  open  literature.  Exam¬ 
ples  of  detailed  studies  are  given  by  Bird  et  al.  [10],  Froment 
and  Bishoff  [16],  Bear  and  Buchlin  [17],  and  Bear  [18]. 

The  mass  transport  equation  needs  to  take  into  account  the 
effect  of  the  medium  porosity.  Moreover,  the  molar  flux  is  mainly 
due  to  convection  in  the  gas  channel,  whereas  in  the  porous 
media,  diffusion  is  the  main  cause.  The  general  form  for  the 
mass  transport  equation  is  [19]: 


g  KXjP) 

RT  dt 


—  —S/Ni  +  (Di 


(24) 


where  s  represents  the  porosity  of  the  electrode,  Xi  the  mole 
fraction  and  Ni  the  mole  flux  of  the  ith  species.  The  reaction  rate 
is  given  by  the  following  expression: 


00; 


—  in  the  TPB 
2  F 

0  elsewhere 


(25) 


In  expression  (25),  F  is  Faraday  constant,  and  j  is  defined 
in  expression  (18).  As  for  the  electrical  problem,  if  the  reaction 
zone  is  reduced  to  the  electrode-electrolyte  interface,  the  reac¬ 
tion  rate  is  zero  in  all  the  electrode,  and  the  mass  flux  is  given 
as  a  boundary  condition. 

A  well  established  way  of  describing  the  diffusion  phe¬ 
nomenon  in  the  SOFC  electrodes  is  through  either  the  Fick’s 
law  [20-24],  or  the  Stefan-Maxwell  equation  [25-28].  Some 
authors  use  more  complex  models,  like  for  example  the  dusty- 
gas  model  [29]  or  derived  from  it  [30,3 1  ] .  A  comparison  between 
the  three  approaches  is  reported  by  Suwanwarangkul  et  al.  [19], 
who  concluded  that  the  choice  of  the  most  appropriate  model 
is  very  case-sensitive,  and  should  be  selected,  according  to  the 
specific  case  under  study. 

If  Fick’s  law  is  considered,  the  mass  flux  is  given  by  the 
following  expression  [19]: 


A,e ff  d(XjP) 
RT  dz 


(26) 


When  the  Stefan-Maxwell  expression  is  used,  the  diffusion 
of  each  gas  species  is  linked  to  the  other  species  composing  the 
gas  mixture  [19]: 


j=hjfi 


XjNj  -  XjNj 

Aj,eff 


J^dXi 
RT  dz 


(29) 


with  n  indicating  the  number  of  gas  mixture  components,  and 
Dij  Q ff  the  effective  binary  diffusion  coefficient. 

A  more  elaborated  expression,  that  takes  into  account  the 
Knudsen  diffusion,  is  given  by  the  dusty-gas  model: 


Ni 


n 


Dik,e  ff 


+ 


E 


XjNj  -  XjNj 
Aj,eff 


P_  fdXA 

RT  V  dz  ) 


(30) 


where  D^,e ff  indicates  the  effective  Knudsen  diffusion  coeffi¬ 
cient  (Veldsink  et  al.  [29]). 

The  momentum  equation  inside  the  porous  media  is  provided 
by  Darcy’s  law: 

k  - 

u  =  — VP  (31) 

T] 

where  k  is  the  permeability  of  the  porous  media,  and  r)  is  the 
dynamic  viscosity. 


2.3.  Current  collectors 

Phenomena  acting  in  the  current  collectors  are  intrinsically 
easier  than  those  occurring  in  the  other  parts  of  the  fuel  cell. 
In  this  component,  in  fact,  only  electrons  flow,  thus  requiring 
only  the  equations  needed  for  modeling  the  electrical  problem 
and  the  temperature  distribution.  Considering  that  the  current 
collectors  are  usually  made  of  very  high  conductive  materials, 
the  relative  ohmic  resistance,  depending  on  the  model  purpose, 
is  usually  ignored.  As  a  consequence,  even  the  heat  generation, 
that  is  due  only  to  the  Joule  effect,  can  be  ignored. 

Regarding  the  electrical  problem,  Eq.  (20)  is  simplified  for 
the  current  collector  as: 

v20  =  0  (32) 


The  current  distribution  can  be  obtained  by: 


where  z  is  the  direction  of  the  electrode  thickness. 


V7  =  0 


(33) 
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(i.e.  the  charge  conservation  inside  the  conducting  media),  cou¬ 
pled  with  Eq.  (16). 

For  the  energy  balance,  Eq.  (21)  is  still  valid,  but  the  source 
term  is  due  only  to  ohmic  losses: 

Sq  =  aV0  •  V0  (34) 

The  current  collector  is  in  direct  contact  with  one  of  the 
electrodes,  with  the  gas,  and  eventually  with  the  surroundings 
(depending  on  the  geometry,  and,  in  the  case  of  the  planar  con¬ 
figuration,  if  this  is  the  bipolar  plate  on  the  top,  on  the  bottom,  or 
inside  the  stack).  In  the  first  case,  heat  is  exchanged  for  conduc¬ 
tion,  while  in  the  second  and  the  third  for  convection.  Finally, 
due  to  the  high  temperature,  radiation  also  plays  an  important 
role  [32-34]. 

2.4.  Electrolyte 

As  for  the  current  collector,  there  is  no  chemical  reaction 
occurring  in  the  electrolyte.  Ions  flow  from  one  side  to  the  other, 
thus  generating  ohmic  losses  and  heat  due  to  the  Joule  effect. 
The  same  equations  used  for  the  current  collector  can  be  used  to 
describe  the  electrolyte  phenomena.  In  this  case,  however,  heat 
is  exchanged  only  for  conduction.  The  equations  to  be  used  are 
here  reported  for  the  reader’s  convenience: 


v20  =  0 

(35) 

< 

II 

O 

(36) 

J  =  — crV0 

(37) 

dQ 

~+VQ  +  Sq=0 

(38) 

Sq  =  <tV0  •  V0 

(39) 

3.  Numerical  model 

3.1.  Equations  discretization 

The  governing  equations  that  model  an  SOFC  are  highly 
non-linear  and  self-coupled,  which  make  it  impossible  to  obtain 
analytically  exact  solutions.  Therefore,  the  equations  must  be 
solved  by  discretization  thus  converting  them  to  a  set  of  numer¬ 
ically  solvable  algebraic  equations. 

The  appropriate  solution  algorithm  to  solve  a  system  of  partial 
differential  equations  strongly  depends  on  the  presence  of  each 
term  and  their  combinations. 

The  three  most  commonly  used  computational  techniques  to 
discretize  a  system  of  partial  differential  equations  are: 

-  Finite  Difference  Method  (FDM)  [35]. 

-  Finite  Volumes  Method  (FVM)  [36,37]. 

-  Finite  Elements  Method  (FEM)  [38-40]. 

Each  approach  allows  transforming  the  continuous  problem 
into  a  discrete  problem  with  a  finite  number  of  degrees  of  free¬ 
dom. 


The  basic  idea  of  the  FDM  is  to  replace  the  individual  terms 
in  partial  differential  equations  with  finite  difference  forms  at  a 
discrete  set  of  points  (the  mesh). 

In  both  the  FEM  and  the  FVM  the  spatial  domain  is  divided 
into  a  finite  number  of  sub-domains  (elements  or  volumes)  and 
the  attention  is  focused  on  these  sub-domains  rather  than  on  the 
nodes  (as  it  is  for  the  FDM).  In  the  FEM  and  the  FVM,  in  fact,  the 
integral  form  of  the  conservation  equations  is  taken  into  account 
and  the  integration  is  made  over  the  prescribed  finite  volumes. 
It  is  important  to  underline  that  the  finite  volume  over  which  the 
integration  is  made  may  coincide  with  a  sub-domain  of  the  mesh. 
Central  to  FVMs  and  FEMs  is  the  concept  of  local  conservation 
of  the  numerical  fluxes.  The  integration  of  the  convective  term  of 
a  conservation  equation  over  a  finite  volume  leads,  in  fact,  to  the 
flux  of  the  generic  variable  through  the  finite  volume  boundary 
surface,  by  applying  the  well-known  Gauss  theorem. 

The  main  difference  between  a  FEM  and  a  FVM  is  that  in 
the  first,  a  function  is  assumed  for  the  variation  of  each  variable 
inside  each  element  while  in  the  latter  this  function  is  always 
equal  to  1 . 

Another  important  feature  of  FVMs  and  FEMs  is  that  they 
may  be  used  on  arbitrary  geometries,  using  structured  or  unstruc¬ 
tured  meshes. 

As  shown  above,  independently  from  the  numerical  method, 
the  spatial  continuum  is  divided  into  a  finite  number  of  discrete 
cells,  and  finite  time-steps  are  used  for  dynamic  problems. 

In  order  to  perform  space  discretization,  the  domain  over 
which  the  governing  equations  apply  is  filled  with  a  predeter¬ 
mined  mesh  or  grid.  The  mesh  is  made  up  of  nodes  (i.e.  grid 
points)  and/or  elements  at  which  the  physical  quantities  (i.e. 
unknowns)  are  evaluated.  Neighboring  points  are  used  to  calcu¬ 
late  derivatives. 

Mesh  generation  is  a  very  complex  task  for  applied  problems 
and  many  different  approaches  to  it  have  been  developed  and 
are  currently  under  study  [41-44]. 

First  of  all,  computational  grids  are  classified  as  structured 
or  unstructured  meshes,  even  if  each  of  these  classes  comprises 
a  broad  list  of  meshing  techniques. 

The  simplest  structured  mesh  is  the  Cartesian  mesh,  where 
nodes  are  distributed  regularly  at  equal  distance  from  one 
another.  Usually  FDMs  are  applied  to  structured  Cartesian 
meshes.  More  generally,  a  grid  is  classified  as  structured  when 
only  the  physical  location  (coordinates)  of  each  node  must  be 
stored  since  the  identity  of  neighboring  mesh  points  is  known 
implicitly.  This  comprises  also  multiblock  structured  grid  gen¬ 
eration  schemes  (a  collection  of  several  structured  blocks  con¬ 
nected  together).  The  main  advantage  of  structured  grids  is 
the  simplicity  in  terms  of  application,  development,  computa¬ 
tion  and  visualization.  The  automatic  connectivity  information 
implies  that  structured  grids  require  the  lowest  amount  of  mem¬ 
ory  for  a  given  mesh  size  and  the  related  simulation  is  faster. 
However,  structured  grids  may  represent  a  severe  limitation  for 
practical  engineering  purposes  especially  when  the  geometry  is 
complex  and  there  is  a  need  for  high  resolution  in  particular 
regions  of  the  domain. 

In  an  unstructured  mesh  each  node  can  have  a  different  num¬ 
ber  of  neighbors  and  elements  have  different  shapes  and  sizes. 
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Therefore,  connectivity  information  must  be  explicitly  defined 
and  stored.  The  unstructured  grid  approach,  that  has  gained  pop¬ 
ularity  with  the  enormous  advancements  of  computer  technol¬ 
ogy,  allows  handling  complex  geometries  with  a  lower  number 
of  elements  and  a  much  easier  realisation  of  local  and  adaptive 
grid  refinement  [44] . 

A  detailed  description  of  mesh  generation  theory  goes  beyond 
the  present  paper  and  the  reader  can  refer  to  literature  [41-44]. 
However,  it  is  important  to  underline  that  the  accuracy  and  the 
stability  of  any  numerical  computation  are  significantly  influ¬ 
enced  by  the  particular  meshing  strategy. 

In  the  case  of  dynamic  (unsteady)  problems,  even  after  the 
space  discretization,  we  still  have  to  solve  a  set  of  ordinary 
differential  equations  in  time.  Therefore,  the  second  step  is  to 
discretize  the  temporal  continuum.  This  is  usually  done  by  a 
finite  difference  approximation  with  the  same  properties  of  an 
FDM  in  space.  Depending  on  the  instant  in  which  the  informa¬ 
tion  is  taken,  the  time-discretization  leads  to: 

-  an  explicit  scheme,  when  the  solution  at  time  step  n  +  1 
depends  only  on  the  known  solution  at  time  step  n  (time¬ 
marching  solution); 

-  an  implicit  scheme,  when  the  time-discretization  step  leads  to 
a  non-linear  algebraic  system  of  equations  that  must  be  solved 
to  calculate  the  solution  at  time  n  +  1 . 

Associated  with  numerical  problems  is  the  concept  of  sta¬ 
bility.  A  numerical  scheme  is  stable  when  a  solution  is  reached 
even  with  large  time-steps  (unsteady  problems)  or  iteration  steps 
(algebraic  system  of  equations  iteratively  solved).  Therefore,  the 
size  of  the  time-step  or  of  the  iteration- step  is  dictated  by  sta¬ 
bility  requirements.  It  must  be  kept  in  mind  that  stability  does 
not  mean  accuracy:  an  implicit  scheme  of  a  dynamic  problem 
is  unconditionally  stable  but  the  solution  obtained  with  large 
values  of  the  time  step  may  not  be  realistic. 

3.2.  Approaches  for  the  problem  solution 

According  to  the  specific  needs  of  the  model,  the  equations 
previously  defined  can  be  implemented,  formulating  some  sim¬ 


plifications  in  the  equations  themselves  or  neglecting  one  or 
more  geometry  dimensions.  In  the  following,  the  formulations 
of  the  problem,  according  to  different  approaches  are  briefly 
presented. 

3.2.1.  Three-dimensional  approach 

This  approach  is  used  when  detailed  information  on  the  fuel 
cell  operation  is  needed.  In  this  case,  all  the  equations  given  in  the 
previous  sections  are  solved.  If  the  dynamic  behavior  of  the  fuel 
cell  is  not  relevant,  the  time-dependent  terms  are  neglected,  thus 
the  results  are  relative  to  the  steady  state  condition.  The  boundary 
conditions  described  in  Section  4  need  to  be  implemented. 

3.2.2.  Two-dimensional  approach 

In  this  case,  one  geometrical  dimension  is  neglected  and  a 
2D  section  is  considered  as  representative  for  the  entire  cell 
operation.  Fig.  4  shows  how  a  micro-tubular  SOFC  geometry  can 
be  reduced  in  a  2D  domain.  Thanks  to  the  perfect  cell  symmetry, 
the  only  assumption  in  this  case  is  that  for  each  angular  section, 
the  cell  behavior  is  the  same. 

For  other  geometries,  a  2D  approach  prompts  to  some 
assumptions  and  simplifications,  causing  a  reduction  in  the 
resulting  information.  Consider,  for  example,  Fig.  5,  where  a 
planar  SOFC  is  depicted.  The  representative  two-dimensional 
section  can  be  chosen  in  several  ways.  In  Fig.  5,  three  cases  are 
identified.  In  case  1,  one  section  along  the  cell  is  chosen.  In  this 
case,  the  section  cannot  be  considered  as  representative  for  the 
entire  cell  operation.  In  fact,  the  gas  species  are  consumed  and 
produced  while  flowing  into  the  channel.  As  a  result,  the  current 
density,  the  temperature  and  all  the  other  physical  proprieties 
vary  in  the  v  direction.  Sometimes  a  quasi-three-dimensional 
approach  is  used  by  discretizing  the  cell  into  a  finite  number  of 
sections  modeled  with  a  two-dimensional  approach  and  dynam¬ 
ically  coupling  (input  of  ith  section  =  output  of  (i  —  l)th  section) 
the  resulting  two-dimensional  models. 

Another  approach  is  to  choose  a  section  in  the  vz  plane  (cases 
2  and  3),  thus  the  variation  in  the  v  direction  can  be  considered. 
However,  in  case  2,  the  section  does  not  comprise  the  gas  chan¬ 
nel,  while,  in  case  3,  the  current  collector  is  not  in  direct  contact 
with  the  electrode.  Therefore,  the  electrical  potential  of  the  exter- 
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Fig.  4.  Simplification  of  a  micro-tubular  geometry  in  2D. 
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Case  1 


Fig.  5.  Possible  cross-sections  for  a  2D  representation. 


nal  electrode  boundary  is  not  a  constant  (contrarily  to  case  2). 
Moreover,  in  both  cases  2  and  3,  the  ohmic  resistance  variation 
along  the  neglected  direction  needs  to  be  taken  into  account  in 
the  model  formulation.  In  fact,  the  general  expression  to  evaluate 
the  ohmic  resistance  is: 

^ohm  =  c  (40) 

a  S 

where  t  and  S  are,  respectively,  the  thickness  and  the  section 
of  the  conducting  medium,  where  the  current  passes  through. 
Expression  (40)  shows  that  the  ohmic  resistance  is  strictly  linked 
to  the  conductive  medium  geometry  and  the  current  path  (i.e.  the 
terms  “thickness”  and  “section”).  For  this  reason,  equivalent  cir¬ 
cuit  models  for  evaluating  the  equivalent  resistance  domain  are 
implemented  (e.g.  [45]).  The  most  complex  situation  is  that  of 
the  Siemens- Westinghouse  design  (Fig.  6).  In  this  configura¬ 
tion,  the  internal  part  of  the  tube  constitutes  the  cathode,  while 
the  external  is  the  anode.  Single  cells  are  connected  to  each  other 
by  conductive  strips  inserted  to  the  anode  of  one  cell  and  touch¬ 
ing  the  cathode  of  the  contiguous  single  cell.  The  current  path  of 
the  cell  is  shown  in  the  cross-section  of  the  cell.  Clearly,  it  is  not 
easy  to  determine  what  t  and  S  are,  according  to  expression  (40), 
thus  requiring  equivalent  circuit  models.  Bharadwaj  et  al.  [46] 
developed  an  analytical  model  to  determine  the  ohmic  polariza¬ 
tion  of  the  Siemems- Westinghouse  tubular  cells.  Using  a  similar 
approach,  in  other  publications,  the  authors  also  computed  the 
ohmic  resistance  of  the  novel  high  power  density  (HPD)  design 
of  Siemens- Westinghouse  [47]. 

Another  possibility,  not  shown  in  Fig.  5,  is  to  choose  the 
x-y  plane  (e.g.  [48]).  However,  in  this  case,  there  is  only  one 
plane  (i.e.  the  x-y  itself)  representing  the  entire  cell,  thus  no 
information  on  single  cell  components  (electrodes,  electrolyte, 
current  collectors)  can  be  directly  deduced. 

3.2.3.  One -dimensional  approach 

In  ID  models  only  one  geometrical  dimension  is  considered 
and  the  fuel  cell  is  represented  by  a  line.  This  is  equivalent 
to  make  the  assumption  that  the  variation  of  the  fluid  and  the 
electrical  properties  along  the  other  two  directions  is  negligible. 
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As  it  is  for  2D  models,  reducing  a  fuel  cell  into  a  ID  domain 
is  a  non-trivial  task  and  the  resulting  simplification  leads  to 
neglecting  important  phenomena  that  must  be  considered  by 
other  means  (i.e.  concentration  losses). 

In  ID  planar  SOFC  models  the  dimension  considered  is 
usually  determined  by  the  gas  flow  direction.  Therefore,  the 
coordinate  varies  along  a  direction  parallel  to  the  axes  of  the 
gas  channels  (referring  to  Fig.  5,  v-direction),  and  this  means 
that  cross-flow  planar  SOFCs  cannot  be  simulated  by  means 
of  ID  models.  By  making  more  simplifying  assumptions  Stan- 
daert  et  al.  [49]  have  proposed  an  analytical  ID  model  for  planar 
molten  carbonate  fuel  cells,  later  extended  to  SOFCs  by  Bove 
et  al.  [50]. 

For  tubular  SOFCs  the  kept  dimension  is  usually  the  tube 
axis  (which  coincides  with  the  direction  of  the  fuel  and  oxidant 
flow),  as  in  the  models  presented  by  Haynes  and  Wepfer  [51], 
by  Ota  et  al.  [52]  and  by  Bove  et  al.  [50],  as  well  as  by  Fi  and 
Suzuki  [53]  and  by  Yokoo  and  Take  [54]  that  extended  the  ID 
analysis  to  the  pre-reformer  and  the  internal  reformer  following 
the  flow  direction. 

Gas  composition  and  flow  rate  in  each  gas  channel  are  average 
values  of  the  channels  section.  In  the  majority  of  the  ID  mod¬ 
els  presented  in  literature  a  simplified  electrochemical  model  is 
considered  and  the  fuel  cell  voltage  is  calculated  by  subtracting 
losses  by  the  ideal  (Nernst)  potential.  The  porosity  effects  of  the 
electrodes  (neglected  in  a  ID  model)  are  considered  by  means 
of  concentration  losses,  thus  Eqs.  (24)-(31)  are  replaced  with  a 
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voltage  reduction  term  (concentration  loss),  that  can  be  modeled 
with  easy  or  more  sophisticated  approaches. 

An  equivalent  electrical  circuit  must  be  constructed  in  order 
to  take  into  account  all  the  ohmic  resistances  (also  those  along 
the  neglected  dimensions)  along  the  current  path. 

This  kind  of  models  allows  to  calculate  current  density  and 
voltage  distribution  along  the  fuel  cell  length.  If  heat  transfer  is 
taken  into  account  temperature  variation  along  the  fuel  cell  can 
be  calculated. 

3.2.4.  Zero -dimensional  approach 

Zero-dimensional  models,  also  called  box  models,  are  the 
simplest  ones.  The  fuel  cell  is  represented  by  only  one  box  and 
spatial  averaging  all  dimensions  is  performed.  Thus,  spatial  vari¬ 
ations  are  not  taken  into  account  and  one  or  more  transformations 
(i.e.  global  mass  and  energy  balances)  are  considered  to  define 
output  variables  from  input  ones.  In  a  dynamic  simulation,  time 
is  the  only  independent  variable. 

A  zero-dimensional  model  could  be  developed  to  examine 
the  impact  of  inlet  composition,  utilization  factors  and  overpo¬ 
tential  on  the  performance  of  an  SOFC  in  terms  of  efficiency  and 
characteristic  curve.  The  thermodynamic  variables,  temperature 
and  pressure,  the  mass  flow  rates  and  the  gas  composition  are 
assumed  to  be  global  values  that  may  change  in  time  (dynamic 
models)  or  from  input  to  output  values. 

However,  considering  that  when  the  spatial  variation  of  the 
variables  is  neglected  (i.e.  the  geometry  does  not  affect  directly 
the  performance  of  the  fuel  cell)  these  models  are  not  suitable  to 
performance  prediction.  More  appropriately,  zero-dimensional 
models  should  be  used  where  attention  is  not  focused  on  the  fuel 
cell  itself  but  on  how  the  fuel  cell  affects  the  performances  of 
the  whole  system.  Two  categories  of  zero-dimensional  models 
are  the  empirical  models  (the  performances  of  the  fuel  cell  are 


already  known  experimentally)  or  the  so-called  “state-of-the- 
art”  components  in  more  complex  systems  (i.e.  combined  energy 
systems). 

Box  models  are  usually  employed,  as  thermodynamic  mod¬ 
els,  for  the  numerical  analysis  of  fuel  cell  based  energy  systems 
(i.e.  SOFC/gas  turbine  hybrid  systems  and  CHP  configurations) 
where  the  single  elements  such  as  compressors,  heat  exchangers, 
fuel  reformer,  partial  oxidizers,  contaminant  removal  appara¬ 
tus,  are  simulated  through  independent  box  models  [55,56].  The 
results  of  each  block  are  the  input  data  for  the  next  block. 

For  instance  Lunghi  and  Ubertini  [57]  integrated  a  zero¬ 
dimensional  model  of  a  planar  SOFC  into  a  system  sequential 
software  (based  on  Aspen  Plus™)  in  order  to  study  different 
configurations  of  SOFC/GT  hybrid  systems  (Fig.  7  shows  one 
power  plant  layout  studied  by  Lunghi  and  Ubertini  [57]).  The 
input  data  to  the  fuel  cell  model  are  gas  temperature,  pressure, 
mass  flow  rate  and  gas  composition.  SOFC  performances  are 
evaluated  through  a  simplified  electrochemistry  model: 


V  —  Vrev  7? ohm  ^act,a  7?act,c  Ocon 


(41) 


where  the  ohmic  resistance  is  estimated  considering  the  resis¬ 
tivity  values  and  the  geometry  of  the  cell  components  and  inter¬ 
connect  materials. 

The  overpotentials  due  to  the  activation  barriers  at  the  elec¬ 
trodes,  77 act,  are  calculated  through  the  Butler- Volmer  Eq.  (18) 
by  means  of  a  mean  current  density  value. 

The  concentration  polarization  losses  are  expressed  as  a  func¬ 
tion  of  the  limiting  current,  i\,  which  is  usually  taken  as  a  measure 
of  the  maximum  rate  at  which  a  reactant  can  be  supplied  to  an 
electrode: 
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Fig.  7.  Low  pressure  SOFC-GT  hybrid  plant  with  steam  production  [57]. 
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Fig.  8.  Variation  of  the  simulated  polarization  curves  if  inlet,  outlet  or  an  average 
chemical  composition  of  the  gases  are  considered  in  a  zero-dimensional  model 
[58]. 

Finally,  the  thermal  behavior  of  the  cell  is  evaluated  through 
the  energy  balance  of  the  system  assuming  one  representative 
temperature  for  each  cell  element  and  reactant  gas. 

The  outputs  of  this  box  model  are  the  fuel  cell  power  output 
and  the  outlet  gas  composition  and  temperature. 

In  zero-dimensional  models  it  is  assumed  that  the  fuel  cell 
is  a  zero-dimensional  space  (a  point).  However,  two  conditions, 
inlet  and  outlet,  are  present  and  an  ambiguity  arises  when  cell 
performances  are  calculated.  Cell  voltage,  in  fact,  depends  on  gas 
concentration  and  different  results  are  obtained  if  inlet  or  outlet 
gas  composition  are  used.  Fig.  8  shows  how  the  characteristic 
curve  of  an  SOFC  changes  using  inlet,  outlet  or  an  average  gas 
concentrations  [58]. 

Fuel  cell  zero-dimensional  models,  as  the  one  described 
above,  are  usually  based  on  assumptions,  parameters  and  practi¬ 
cal  information  provided  in  literature  or  taken  from  experiments. 
The  objectives  of  this  kind  of  analyses  could  be  the  evalua¬ 
tion  of  optimal  power  plant  design,  configuration,  lay-out  and 
components  sizing  and  of  the  related  performances  in  terms  of 
efficiency  and  power  output.  Parametric  studies  of  SOFC  based 


energy  systems  can  be  performed  in  a  relatively  short  time  for 
optimizing  the  performances  and  the  operation  in  function  of 
design  and  operation  parameters  of  the  fuel  cell  (i.e.  stack  size, 
utilization  factors,  cell  temperature)  and  of  the  other  system 
components. 

For  instance,  Costamagna  et  al.  [55]  developed  a  zero¬ 
dimensional  model  of  a  SOFC/micro  gas  turbine  (MGT)  hybrid 
system  to  analyse  a  SOFC/MGT  hybrid  system  under  variable 
fuel  cell  and  MGT  operating  conditions,  with  particular  empha¬ 
sis  on  the  MGT  rotational  speed.  The  fuel  cell  group  is  a  tubular 
SOFC  modeled  through  macroscopic  equations  that  express  the 
balance  between  inlet  and  outlet  conditions.  The  system  simula¬ 
tions  have  been  performed  on  the  basis  of  experimental  maps 
for  the  compressor  and  the  turbine  and  of  indications  from 
the  literature  on  SOFC  and  MGT  plants.  This  simplified  zero¬ 
dimensional  approach  was  used  in  order  to  save  computational 
time  while  maintaining  an  acceptable  accuracy  in  the  overall 
plant  performances  evaluation. 

As  a  result  they  could  evaluate  the  performances  of  the 
SOFC/GT  energy  system  under  part-load  and  off-design  con¬ 
ditions,  by  varying  a  huge  quantity  of  parameters.  In  particular, 
Fig.  9  shows  the  hybrid  plant  performances,  in  terms  of  effi¬ 
ciency  versus  net  non-dimensional  power  (part-load  conditions), 
at  different  SOFC  air  utilization  factors  and  current  densities 
under  variable  MGT  speed.  The  results  show  that  the  possibility 
of  varying  the  MGT  rotational  speed  allows  maintaining  high 
overall  system  efficiency  even  at  very  low  load  conditions. 

This  kind  of  analysis  is  of  fundamental  importance  to  estimate 
the  potential  of  the  hybrid  system  for  distributed  energy  that 
is  characterized  by  highly  variable  electric  and  thermal  loads. 
Combining  the  graph  shown  in  Fig.  9  with  the  load  request,  in 
fact,  would  allow  to  find  the  optimal  operating  parameters  for 
the  chosen  application. 

3.2.5.  Multi-dimensional  approach 

It  is  a  common  practice  to  combine  more  approaches,  espe¬ 
cially  in  complex  energy  systems  as  fuel  cells.  For  instance,  a 
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Fig.  9.  Part-load  performances  of  a  SOFC/MGT  hybrid  system.  MGT  stands  for  micro  gas  turbine  and  P/Pd.p  stands  for  net  non-dimensional  power  [55]. 
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zero-dimensional  model  may  be  used  to  estimate  initial  values 
for  a  three-dimensional  simulation.  Moreover  a  one-dimensional 
approach  can  be  used  to  model  the  flow  field  in  the  gas  channels 
together  with  a  two-dimensional  model  of  the  temperature  and 
the  current  paths  at  the  electrodes. 

Multi-dimensional  models  could  be  considered  as  modular 
softwares  where  the  different  components  could  be  simulated 
through  OD,  ID,  2D  or  even  3D  models. 

4.  Boundary  conditions  (internal  and  external) 

The  accuracy  of  a  modeling  technique  prediction  is  highly 
sensitive  to  the  specification  of  the  initial  and  the  boundary 
conditions,  which  usually  represents  the  major  part  of  the  com¬ 
putational  effort.  In  a  few  words,  boundary  conditions  are  the 
mathematical  description  of  the  different  situations  that  pro¬ 
duce  different  results  within  the  same  physical  system  (same 
governing  equations).  A  proper  specification  of  the  boundary 
conditions  is  necessary  to  allow  the  calculation.  Once  boundary 
conditions  are  defined,  the  so-called  “properly-posed  problem” 
is  reached.  Moreover,  it  must  be  noted  that  in  modeling  an  SOFC 
there  are  various  intrinsically  coupled  flow,  thermal,  chemical 
and  electrochemical  phenomena  and  different  physical  systems 
to  be  modeled.  This  means  that  the  definition  of  the  boundary 
conditions  incorporates  both  external  boundary  conditions  and 
interface  boundary  conditions. 

The  type  of  partial  differential  equations  and  the  dis¬ 
cretization  approach  determine  the  form  of  the  boundary  con¬ 
ditions,  that  are  usually  classified  either  in  terms  of  their 
mathematical  description  or  in  terms  of  the  physical  type  of 
the  boundary. 

For  steady  state  problems  the  following  two  types  of  spatial 
boundary  conditions  are  identified: 

-  Dirichlet  boundary  condition,  when  the  generic  variable  on 
the  boundary  assumes  a  known  and  constant  value; 

-  Neuman  boundary  condition,  when  the  derivatives  of  the 
generic  variable  on  the  boundary  are  known  and  this  yields 
an  extra  equation. 

Physical  boundary  conditions  depend  on  the  specific  problem 
to  be  solved. 

In  the  fluid  problem  within  the  gas  channels  the  following 
boundary  conditions  are  observed: 

-  stationary  solid  walls  (sidewalls)  where  the  velocity  compo¬ 
nents  vanish  at  the  walls  and  the  heat  flux  must  be  defined; 

-  interface  between  the  gas  channel  and  the  porous  media  where 
mass  and  energy  crosses  the  boundary  surface  in  either  direc¬ 
tion;  these  conditions  are  therefore  usually  given  imposing  the 
continuity  of  mass  and  energy  fluxes; 

-  open  flows  (inlet  and  outlet)  where  the  fluid  enters  or  leaves 
the  domain;  fluid  velocity  or  pressure  values  together  with,  or 
the  mass  flow  rate  should  be  defined  as  well  as  the  specific 
enthalpy  for  the  energy  equation;  the  initial  molar  or  mass 
fraction  of  the  different  species  must  be  also  imposed  at  inlet. 


In  the  following  the  possible  boundary  conditions,  for  each 
previously  defined  equation,  are  given. 

4.1.  Channels  flow 

In  order  to  solve  Eqs.  (6),  (10)  and  (12),  a  number  of  uncou¬ 
pled  conditions  equal  to  the  number  of  unknowns  need  to  be 
specified  at  the  boundaries  (boundary  conditions).  These  are  the 
mass  fractions  T/  (continuity  equation),  the  fluid  temperature  T 
(energy  equation)  and  a  characteristic  of  the  flow  field  (momen¬ 
tum  equation),  i.e.  u  or  P.  Taking  as  a  reference  the  case  of  Fig.  1 
(planar  configuration),  the  molar  fraction  needs  to  be  specified 
at  the  inlet,  while  at  the  channel  boundaries,  there  is  no  varia¬ 
tion.  At  the  electrode/channel  flow  interface  and  at  the  outlet, 
the  continuity  is  imposed.  The  velocity  of  the  fluid  is  zero  at  the 
channel  walls,  while  it  is  specified  at  the  inlet  and  outlet.  Alter¬ 
natively,  it  can  be  convenient  to  specify  (usually  at  the  outlet) 
the  pressure. 

4.2.  Electrodes 

According  to  the  schematization  of  Fig.  1,  the  electrodes  are 
in  contact  with  the  electrolyte  on  one  side,  and  with  the  gas  and 
the  current  collector  on  the  other  side.  Although  the  tubular  con¬ 
figuration  presents  a  different  geometry,  it  is  always  possible  to 
identify  these  three  boundary  faces  [15].  For  the  sake  of  sim¬ 
plicity,  the  planar  configuration  is  taken  as  reference,  and  the 
electrode  domain  is  reported  in  Fig.  10;  however,  in  complete 
analogy,  similar  boundary  conditions  can  be  provided  for  the 
tubular  or  monolithic  configuration. 

Surfaces  SI,  S2,  S4  and  the  other  perimeteral  surface  not 
represented  in  Fig.  10  (because  hidden)  are  insulated,  thus  the 
boundary  conditions  for  the  electrical,  thermal  and  mass  trans¬ 
port  problems  are: 

j'H  =  0  (43) 

Q-h  =  -AVT  •  h  =  0  (44) 

Vm/  •  h  =  0  (45) 

Alternatively,  an  estimation  of  the  heat,  mass  and  current  loss 
can  be  provided.  The  estimation  of  these  quantities  are  usually 
affected  by  a  large  uncertainty  degree,  thus  the  ideal  case  is 
usually  considered. 


S2  - S3 

Fig.  10.  Schematic  representation  of  the  electrode. 
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Surface  S3  represents  the  boundary  between  the  electrode  and 
the  electrolyte.  As  mentioned  in  Section  2.2,  since  the  TPB  is 
usually  very  thin  and  close  to  the  electrode-interface  boundary, 
the  TPB  can  be  assimilated  with  surface  S3.  Eqs.  (17),  (20),  (22) 
and  (25)  reduce,  respectively  to: 


V7  =  0  (46) 

V20  =  0  (47) 

Sq  =  ctV0  •  V0  (48) 

on  =  0  (49) 


In  this  case,  the  electrochemical  reactions  (3)  and  (4)  are 
considered  to  occur  at  boundary  S3,  thus  the  boundary  condition 
for  the  electrical  problem  is  provided  by  Eq.  (18)  that  are  here 
reported  for  the  reader’s  convenience: 


J  -n  =  jo 


exp 


Oi]Jhct\ 

RT  ) 


exp 


\  RT  ) 


The  condition  for  the  thermal  problem  is: 

Q  -  n  =  ——  T A  S  +  Jr/act 
2  F 


while,  for  the  mass  transport,  the  conditions  are: 
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In  expressions  (52)-(54),  the  symbol  M/  represents  the 
molecular  weight  of  the  ith  species.  Expressions  (52)  and  (53) 
are  referred  to  the  anode,  while  expression  (54)  to  the  cathode. 

Surface  S5  is  the  boundary  between  the  current  collector  and 
the  electrode.  The  following  conditions  are  applied: 

J  •  n  =  (1  •  n).  (55) 

V  /  interconnect  v  7 

Q  •  n  =  (q  -  it)  (56) 

V  /  interconnect 

rhi  •  n  =  0  (57) 


The  first  two  expressions  state  the  current  and  heat  flux  conti¬ 
nuity,  while  the  third  states  the  mass  insulation,  due  to  the  current 
collector. 

Finally,  surface  S6  is  in  direct  contact  with  the  gas.  In  this 
case  there  is  a  mass  continuity  with  the  channel  flow,  insulation 
for  the  electric  current,  and  a  heat  flux  due  to  convection  with 
the  gas  and  radiation  from  the  gas  channel  walls: 

J  -n  =  0  (58) 

Q  '  Tl  =  /z(7s6  —  Tgas)  +  £2 rad-surf  (59) 

(fHi  ’  ^Ochannel  flow  (60) 

In  expression  (59),  h  is  the  convective  heat  transfer  coeffi¬ 
cient  between  the  gas  flow  and  the  electrode,  and  Grad-surf  is  the 


S5 

Fig.  11.  Current  collector/gas  channel  of  a  flat  planar  SOFC. 


radiation  heat  transferred  between  S6  and  any  other  visible  sur¬ 
face.  Murthy  and  Fedorov  [33]  noted  that  the  surface-to-surface 
approach  (as  that  of  relation  (59))  could  lead  to  some  tempera¬ 
ture  prediction  mistakes.  More  accurate  results  can  be  expected 
considering  the  absorption,  emission  or  scattering  in  the  media. 

4.3.  Current  collectors 

Fig.  11  schematically  represents  the  current  collector  of  a 
planar  SOFC.  The  surfaces  that  are  thermally  and  electrically 
insulated  are  S 1 ,  S2,  S3  and  forth  perimeteral  face  (not  indicated 
Fig.  11).  For  these  faces,  the  following  conditions  apply: 

J-n  =  0  (61) 

Q  n  =  -kVT  •  n  =  0  (62) 

Surface  S5  is  electrically  insulated,  thus  expression  (61) 
applies,  but  convective  and  radiation  heat  transfer  need  to  be 
considered. 


Q  •  n  —  h(T$ 5  Tgas)  T  Grad-surf  (63) 

For  surface  S4,  the  continuity  conditions  needs  to  be  applied: 
j  -  n  =  (j-n).  (64) 

V  /  interconnect  v  7 

Q  ■  n  =  (Q  -  n)  (65) 

\  /  electrode 

Similarly  to  relation  (64)  and  (65)  the  continuity  conditions 
for  boundary  S6  (hidden  in  Fig.  11),  are  applied. 


4.4.  Electrolyte 


Referring  to  Fig.  1,  the  electrolyte  can  be  regarded  as  a  par¬ 
allelepiped,  where  the  surfaces  constituting  the  perimeter  are 
insulated: 


./•/?=  0  (66) 

Q  n  =  — AVT  •  n  =  0  (67) 
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For  the  surfaces  at  the  top  and  the  bottom,  continuity  with 
the  two  electrodes  needs  to  be  expressed. 


Q  h  =  ( Q  •  h) 

V  /  electrode 

J  •  n  =  0.5(7  •  n)  ,  „  , 

V  /  electrode 


(68) 

(69) 


The  term  0.5  in  Eq.  (69)  takes  into  account  that,  as  stated  in 
Eq.  (4),  every  two  electrons  released  (reacting)  at  the  anode 
(cathode),  one  02~  ion  flows  through  the  electrolyte.  As  a 
reminder,  Eqs.  (35)-(37)  assume  that  the  ionic  current  follows 
the  ohmic  law. 


5.  Validation 

Since  the  numerical  result  of  a  computation  is  only  an  approx¬ 
imation  of  real  life  and  considering  that  convergence  is  not  a 
sufficient  condition  to  stability,  a  necessary  step  in  the  devel¬ 
opment  of  a  computational  model  is  validation.  Validation  is 
performed  through  the  comparison  between  numerical  results 
and  experimental  data  to  establish  the  range  of  validity  and  the 
accuracy  of  the  code.  Blind  acceptance  of  computed  results  are 
not  a  good  basis  to  make  engineering  decisions  involving  finan¬ 
cial  and  time  schedule  risks.  However,  it  is  important  to  underline 
that  even  a  validation  against  many  and  the  most  complex  prob¬ 
lems  cannot  prove  that  the  model  is  totally  error-free. 

If  poor  agreement  is  found  between  experiments  and  compu¬ 
tations,  the  task  is  to  understand  error  sources  both  in  the  code 
and  in  the  measuring  system:  the  input  data  may  involve  too 
much  guess-work  or  imprecision;  the  available  computer  power 
may  be  too  small  for  high  numerical  accuracy;  the  scientific 
knowledge  base  may  be  inadequate. 

6.  Results 

In  the  present  section,  results  obtained  from  the  approaches 
described  in  Section  3.2  are  provided. 

6.1.  Three-dimensional  approach 

Three-dimensional  models,  coupling  and  solving  electrical, 
fluid  dynamic,  chemical  and  thermal  problems  are  very  com¬ 
plex,  time  consuming  and  can  require  strong  numerical  simula¬ 
tion  expertise  and  hardware.  For  this  reason  3D  models  usually 
found  in  the  literature  model  only  some  of  the  above  problems, 
and  neglect  others.  However,  due  to  the  recent  improvements 
of  the  simulation  tools  (hardware  and  software),  some  authors 
have  implemented  all  the  equations  (even  if  with  some  simpli¬ 
fied  assumptions).  Fig.  12  shows  the  mesh  of  a  five  channels, 
cross-flow,  electrolyte  supported  planar  SOFC  [59].  The  related 
current  density  distribution  is  depicted  in  Fig.  13.  The  details  of 
the  results  and  the  relevant  amount  of  information  that  can  be 
deducted  from  this  analysis  is  clear  from  Fig.  13.  Such  results 
can  be  obtained  for  all  the  other  physical  characteristics  con¬ 
sidered  in  the  model,  such  as  temperature,  velocity  distribution, 
concentration  distribution,  electrical  potential,  etc. 


Anode  + 
Inlets 


Anode  < 
Channels 


Anode 

Outlets 


Cathode 

Channels 


Cathode 

Inlets 


Cathode 

Outlets 


Fig.  12.  Mesh  of  a  five  channels,  cross-flow,  electrolyte-supported  SOFC 
[59]. 


6.2.  Two-dimensional  approach 

As  explained  in  Section  3.2,  there  are  different  approaches  to 
simplify  a  domain  from  3D  to  2D.  In  the  present  section,  as  an 
example,  the  case  of  a  micro-tubular  SOFC  is  illustrated.  Fig.  4 
depicts  the  micro-tubular  SOFC,  modeled  by  Bove  and  Sammes 
[14].  In  this  configuration,  the  current  collectors  are  placed  at 
the  two  ends  of  the  tube,  thus  representing  two  caps  for  the  tube. 
One  cap  touches  the  outer  part  of  the  tube  (cathode)  and  the 
other  the  inner  part  (anode).  Due  to  the  circular  symmetry  of 
the  cell,  modeling  the  fuel  cell  in  a  2D  environment,  does  not 
lead  to  any  loss  of  information.  Results  for  the  current  path  in 
the  cell  are  depicted  in  Fig.  14  [14].  Please  note  that  since  the 
electrolyte  and  cathode  are  much  thinner  than  the  anode,  only  the 
current  path  trough  the  anode  is  visible  in  the  figure.  However, 
since  the  anode  conductivity  is  much  higher  than  that  of  the 
electrolyte,  ions  migration  through  the  electrolyte  is  limited  to 
the  tube  part  very  close  to  the  cathodic  current  collector  (left  hand 
side  of  Fig.  14).  As  a  consequence,  even  electrons  migration 
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Fig.  14.  Current  path  in  a  micro-tubular  SOFC,  with  the  current  collectors  at  the  two  tube’s  ends  [14]. 


along  the  cathode  length  is  negligible.  For  more  details  about  the 
current  path  in  this  configuration,  please  refer  to  [14].  The  model 
allows  to  fully  understand  the  reason  for  a  poor  performance, 
when  current  collectors  are  placed  in  this  configuration.  The 
current  path,  in  fact,  is  mainly  along  the  cell,  thus  a  high  ohmic 
resistance  is  generated.  In  the  same  study,  Bove  and  Sammes  [14] 
also  analyzed  the  performance  achievable  if,  together  with  the 
two  caps,  the  current  is  collected  using  a  wire  wrapped  around 
(internally  and  externally)  the  cell.  The  results  for  the  current 
density  path  is  depicted  in  Fig.  15.  Comparing  Figs.  14  and  15,  it 
is  clear  how  a  reliable  model  can  help  the  research  in  designing 
single  components,  as  well  single  cells. 


6.3.  One-dimensional  approach 

As  fully  explained  in  Section  3.2,  the  use  of  a  one¬ 
dimensional  approach  implies  to  neglect  those  phenomena 
occurring  in  the  two  directions  not  considered  in  the  model. 
More  evolved  models,  however,  even  if  consider  the  gas 
moving  in  one  direction,  take  into  account  the  interac¬ 
tion  between  different  layers  that  constitute  the  single  cell. 
Fig.  16  shows  how  Haynes  and  Wepfer  simplified  a  Siemens- 
Westinghouse  tubular  cells  [60],  and  the  related  results  for  the 
variation  of  the  fluid  temperature  along  the  cell  are  represented 
in  Fig.  17. 


a  Cell  width  (m) 


Fig.  15.  Results  for  current  density  distribution  in  a  micro-tubular  SOFC  with  wires  wrapped  around  the  tube  [14]. 
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Fig.  16.  Schematic  representation  in  ID  of  the  Siemens- Westinghouse  tubular 
cell  [60]. 


Oxidant  Temperature  Profiles 
Air  Supply  Pipe  and  Annulus 


chemical  composition  and  temperature  of  inlet  and  outlet  gas, 
cell  voltage,  power  density  and  efficiency.  Typically,  the  infor¬ 
mation  required  from  these  kind  of  models  are  the  polarization 
curve,  the  power  density  and  the  efficiency  variations  with  the 
current  density  variation.  Fig.  18  represents  the  polarization 
curve  for  a  tubular  Siemens- Westinghouse  SOFC,  operating  at 
different  temperatures. 

7.  Conclusions 

This  paper  focuses  on  the  development  of  a  detailed  numer¬ 
ical  model  of  an  SOFC.  A  complete  analysis  of  the  phenomena 
acting  in  solid  oxide  fuel  cells  is  presented  and  the  three- 
dimensional  mathematical  model  of  each  element  of  the  SOFC, 
built  on  the  basis  of  conservation  and  constitutive  laws,  is 
written.  The  mathematical  model  is  a  complete,  3D  and  time- 
dependent  model  independent  of  the  fuel  cell  geometry  (i.e. 
planar  and  tubular,  monolithic)  and  the  modeling  approach  (i.e. 
time-dependent,  3D,  2D).  Then  the  discretization  to  produce  a 
numerical  analogue  of  the  equations  is  discussed.  The  boundary 
conditions,  that  strongly  influence  the  accuracy  of  the  model  and 
usually  represent  a  significant  part  of  the  computational  effort, 
necessary  to  reach  a  properly-posed  problem  are  shown. 

Finally,  a  literature  review  is  conducted.  Some  results  from 
the  models  previously  developed  by  the  authors  or  found  in  the 
literature  survey  are  presented.  Each  of  them  can  be  reviewed  as 
a  simplified  version  of  the  model  presented  in  the  present  paper. 
The  simplified  assumptions  and  the  related  omissions  that  each 
of  those  approaches  imply  are  also  analyzed. 


Inside  Air  Supply  Pipe 
Inside  Annular  Region 
Inside  ASP  (Bessette  et  al.) 

Inside  Annular  (Bessette  et  al.) 

Fig.  17.  Fluid  temperature  profile  inside  a  Siemens- Westinghouse  tubular  cell 
[60]. 

6.4.  Zero -dimensional  approach 

Results  for  zero-dimensional  models  can  provide  only 
information  on  the  overall  performance,  like  for  example, 


Current  Density  (A/cmA2) 


Fig.  18.  Polarization  curve  of  a  tubular  Siemens- Westinghouse  single  cell  oper¬ 
ating  at  different  temperature,  obtained  through  a  zero-dimensional  model. 
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